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A  one-dimensional,  two-phase,  transient  PEM  fuel  cell  model  including  gas  diffusion  layer,  cathode  cat¬ 
alyst  layer  and  membrane  is  developed.  The  electrode  is  assumed  to  consist  of  a  network  of  dispersed 
Pt/C  forming  spherically  shaped  agglomerated  zones  that  are  filled  with  electrolyte.  Water  is  modeled  in 
all  three  phases:  vapor,  liquid  and  dissolved  in  the  ionomer  to  capture  the  effect  of  dehydration  of  the 
ionomer  as  well  as  flooding  of  the  porous  media.  The  anode  is  modeled  as  a  sophisticated  spatially  reduced 
interface.  Motivated  by  environmental  scanning  electron  microscope  (ESEM)  images  of  contact  angles  for 
microscopic  water  droplets  on  fibers  of  the  gas  diffusion  layer,  we  introduce  the  feature  of  immobile  sat¬ 
uration.  A  step  change  of  the  saturation  between  the  catalyst  layer  and  the  gas  diffusion  layer  is  modeled 
based  on  the  assumption  of  a  continuous  capillary  pressure  at  the  interface.  The  model  is  validated  against 
voltammetry  experiments  under  various  humidification  conditions  which  all  show  hysteresis  effects  in 
the  mass  transport  limited  region.  The  transient  saturation  profiles  clearly  show  that  insufficient  liquid 
water  removal  causes  pore  flooding,  which  is  responsible  for  the  oxygen  mass  transport  limitation  at 
high  current  density  values.  The  simulated  and  measured  current  responses  from  chronoamperometry 
experiments  are  compared  and  analyzed. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

It  is  known  that  the  cathode  of  a  PEM  fuel  cell,  using  hydro¬ 
gen  as  fuel,  is  the  limiting  component  for  high  performance.  The 
slow  kinetics  of  the  oxygen  reduction  reaction  (ORR)  lead  to  a 
strong  decrease  of  the  cell  voltage  in  the  low  current  density  region. 
Additionally,  mass  transport  limitations  and  instability  of  the  cell 
voltage  occur  in  the  high  current  density  region,  where  typically 
liquid  water  blocks  the  pathway  to  the  active  sites.  Therefore,  accu¬ 
mulation  and  transport  of  liquid  water  is  a  major  factor  in  the 
operation  of  a  low-temperature  PEM  fuel  cell. 

Water  is  typically  introduced  through  both  anode  and  cathode 
gas  streams,  additionally  produced  by  the  ORR  at  the  cathode  cat¬ 
alyst  layer  (CCL)  and  forced  to  the  CCL  by  electro-osmotic  drag 
from  the  anode  side.  An  excess  of  liquid  water  in  the  porous  media, 
namely  the  catalyst  layer  and  the  gas  diffusion  layer  (GDL),  blocks 
a  fraction  of  the  open  pores.  This  accumulated  water  reduces  the 
path  available  for  the  transport  of  gaseous  oxygen  or  forces  oxygen 
to  dissolve  and  diffuse  through  water  to  reach  the  active  sites.  This 
phenomenon  is  called  flooding  and  it  is  most  problematic  in  the 
cathode. 
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Dissolved  water  in  the  ionomer  is  essential  for  optimal  proton 
conductivity  in  the  catalyst  layer  (CL)  and  membrane.  The  conduc¬ 
tivity  of  the  ionomer,  which  consists  of  a  fluorocarbon  polymer 
backbone  with  chemically  bonded  sulfonic  acid  groups  as  side 
chains,  is  a  very  strong  function  of  its  water  content  [1  ].  Dehydra¬ 
tion  of  the  ionomer  strongly  increases  the  resistance  to  the  proton 
transport,  causing  ohmic  losses.  Additionally,  the  active  area  in  the 
CL  depends  on  the  amount  of  dissolved  water  in  the  ionomer  due 
to  need  for  a  three-phase  boundary  for  the  reaction.  This  means 
that  the  catalyst  with  the  carbon  support  needs  good  coupling  to 
the  protonic  phase,  which  occurs  only  if  the  water  content  is  cor¬ 
rect.  Therefore,  increasing  the  power  density  of  a  PEM  fuel  cell  is 
basically  a  water  management  problem.  At  present,  many  transport 
phenomena  during  fuel  cell  operation  cannot  be  directly  observed 
or  measured.  Thus,  inverse  modeling  is  a  powerful  tool  to  gain  qual¬ 
itative  insights  into  the  dynamics  of  liquid  water  transport  and  its 
effect  on  fuel  cell  performance. 

Much  model  development  and  simulation  work  has  been  done 
in  recent  years  concerning  water  distribution  in  a  PEM  fuel  cell 
to  analyze  the  corresponding  loss  mechanisms  [2-11].  The  mod¬ 
els  differ  strongly  in  the  complexity  of  the  applied  physics,  the 
examined  components  and  dimensions  of  the  modeling  domains. 
One-dimensional  models  are  widely  used  for  testing  new  model 
approaches  because  the  governing  equations  can  be  implemented 
quickly  and  the  computing  time  required  for  parameter  studies  or 
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time-dependent  problems  is  short  [12-17].  One  part  of  the  water 
management  investigations  are  focused  on  water  transport  in  the 
membrane  forced  by  diffusion,  convection  and  electro-osmotic 
drag.  Investigations  have  been  made  concerning  the  dependence 
of  the  water  content  distribution  on  the  cell  current  or  its  tran¬ 
sient  behavior  while  undergoing  a  load  change  [18].  Wu  et  al.  [17] 
developed  a  ID,  transient  PEM  fuel  cell  model  to  analyze  dynamic 
characteristics  corresponding  to  various  changes  in  working  con¬ 
ditions,  such  as  relative  humidity  or  cell  voltage.  The  phenomena 
of  hydration/dehydration  is  investigated  using  different  membrane 
types.  Overshoots  and  undershoots  are  observed,  whereas  the 
model  does  not  account  for  liquid  water.  Transient  analysis  of  the 
water  content  in  the  membrane  was  also  investigated  by  Yan  et  al. 
[19]  with  an  isothermal,  one-phase  model.  Chen  et  al.  [20]  include 
the  effect  of  membrane  swelling  in  the  transient  analysis  of  the 
water  transport  in  the  membrane  (vapor  and  dissolved)  and  its 
influence  on  performance.  The  influence  of  liquid  water  in  the  CL 
and  GDL  on  the  performance  and  the  interaction  of  liquid  water  in 
the  porous  media  with  the  dissolved  water  in  the  membrane  were 
not  considered  in  these  models.  With  regard  to  liquid  water,  Eiker- 
ling  [21]  developed  a  detailed  cathode  model  focused  only  on  the 
catalyst  layer  (GDL  and  membrane  are  not  considered).  He  analyzed 
the  impact  of  the  electrode  structure  and  its  physical  properties  like 
wettability  on  the  liquid  water  distribution.  In  doing  so,  he  made 
the  assumption  that  the  electrode  consists  of  a  network  of  agglom¬ 
erates,  similar  to  [22-25],  highlighting  the  role  of  a  bimodal  pore 
size  distribution  for  the  water  distribution.  The  relative  proportion 
of  primary  and  secondary  pores  and  their  wetting  properties  deter¬ 
mine  the  interplay  between  reaction  and  transport  processes.  With 
strong  simplifications,  analytical  solutions  are  shown  for  three  dis¬ 
tinct  states:  a  dry  state,  an  optimal  wetting  state  and  a  fully  flooded 
state.  Transient  analysis  was  not  done.  Madhusudana  et  al.  [26] 
present  a  dynamic  agglomerate  model,  but  neglect  the  effect  of 
liquid  water.  Overshoot  and  undershoot  behavior  of  the  current 
density  to  voltage  step  changes  are  attributed  to  species  transporta¬ 
tion  in  the  agglomerate  structure. 

The  influence  of  the  cell  geometry  concerning  e.g.  channel  or 
shoulder  width/length  on  the  saturation  distribution  in  the  porous 
media  was  analyzed  in  two-dimensional  or  three-dimensional  CFD 
models,  among  others,  by  [27,28,3,29,18].  Siegel  et  al.  [29]  used 
an  agglomerate  approach  for  the  catalyst  layer  structure  in  a  2D, 
two-phase  model.  They  investigated  the  saturation  profile  in  the 
porous  media  along  the  channel  and  its  effect  on  the  cell  perfor¬ 
mance.  Therefore,  the  blockage  of  the  active  area  by  liquid  water 
(saturation  s)  is  modeled  by  the  factor  (1  -s)  in  the  expression 
for  the  reaction  kinetics.  Wang  and  Wang  [18]  investigate  the  fuel 
cell  performance  under  load  changes  with  a  3D,  along-the-channel 
model.  The  water  content  distribution  is  analyzed  in  terms  of  the 
voltage  response.  The  model  accounts  also  for  the  catalyst  layer, 
but  neglects  liquid  water  transport  in  the  porous  media.  A  three- 
dimensional  analysis  of  the  distribution  of  reactants  and  products 
under  the  channel  and  particularly  under  the  rib  is  done  by  Berning 
and  Djilali  [3].  For  simplification,  the  catalyst  layer  is  assumed  to 
be  a  spatially  non-resolved  interface  between  the  membrane  and 
GDL.  Maximum  saturation  values  in  the  cathode  GDL  of  about  10% 
at  a  current  density  of  1.2  A  cm-2,  located  under  the  rib  towards  the 
CL,  are  predicted  by  the  model.  The  capillary  pressure-saturation 
relationship  pc[s]  is  modeled  by  the  Leverett  J-function  [30],  which 
is  the  most  common  approach  in  two-phase  fuel  cell  models.  This 
approach  leads  to  a  strong  driving  force  for  the  capillary  transport 
of  liquid  water,  resulting  in  a  low  saturation  level  in  the  case  of  a 
fixed  saturation  of  zero  at  the  boundary  towards  the  gas  channel.  An 
alternative  expression  for  the  capillary  pressure-saturation  rela¬ 
tionship  was  suggested  by  Natarajan  and  Nguyen  [27,28],  extracted 
by  fitting  an  empirical  function  to  the  measured  pc[s]-curve  of  a 


GDL.  The  corresponding  liquid  water  diffusion  obtained  with  this 
expression  is  four  orders  of  magnitude  smaller,  which  means  that 
extremely  high  saturation  gradients  would  be  necessary  in  order  to 
induce  a  flux  of  liquid  water  [3].  Thus,  the  simulation  results  with 
a  3D  cathode  model  show  high  saturation  levels  close  to  1.  A  capil¬ 
lary  pressure-saturation  relationship  was  also  measured  by  Acosta 
et  al.  [31]  and  fitted  to  an  empirical  function.  They  distinguish 
between  the  imbibition  and  drainage  curve,  leading  to  a  satura¬ 
tion  level  of  about  0.11  and  0.5,  respectively,  for  a  conventional 
gas  distributor  modeled  in  2D.  Kumbur  et  al.  [32]  investigated  the 
effectiveness  of  the  Leverett  approach.  They  discuss  the  existence 
of  immobile  saturation,  a  threshold  below  no  continuous  liquid 
phase  is  formed,  and  compare  different  approaches  for  the  rela¬ 
tive  permeability  (Wyllie,  Corey,  Brooks-Corey  and  Van  Genuchten 
model).  They  also  suggest  the  existence  of  a  saturation  discon¬ 
tinuity  between  two  porous  media  with  different  porosities  and 
wettabilities  (GDL/MPL),  discuss  the  characteristics  of  a  mate¬ 
rial  with  mixed  wettability  and  analyze  the  influence  of  contact 
angles  on  the  saturation  distribution.  A  polynomial  fit  of  the  pc[s]~ 
function  was  made  based  on  experimental  pressure  curves.  Nam 
and  Kaviany  [33]  had  also  suggested  immobile  saturation  several 
years  previously.  They  investigated  the  tortuosity  effect  for  the  case 
of  liquid  water  present  in  a  GDL  and  analyzed  the  enhancement 
of  the  saturation  level  with  different  material  layers  of  different 
wettabilities. 

A  comprehensive  one-dimensional  transient  fuel  cell  model 
which  accounts  for  liquid  saturation  in  the  GDL  and  CL,  on  both 
the  anode  and  cathode  sides,  was  developed  by  Ziegler  et  al.  [34]. 
The  model  was  expanded  from  the  steady-state  model  of  Weber  and 
Newman  [35,36]  for  time-dependent  analysis.  Both  models  account 
for  the  Schroeder  paradox,  describing  that  the  membrane  water 
uptake  is  less  when  exposed  to  100%  relative  humidity  (RH)  than 
liquid  water.  Ziegler  et  al.  compared  their  model  against  potential 
sweep  experiments.  Similar  sweep  experiments,  presented  in  [37], 
are  used  to  validate  a  ID  transient  model  developed  by  Shah  et  al. 
[38].  The  model  takes  into  account  the  agglomerate  structure  of 
the  catalyst  layers,  water  in  the  forms  of  vapor,  liquid  and  dissolved 
in  the  ionomer.  A  feature  of  this  model  is  the  boundary  condition 
for  saturation  at  the  interface  GDL  channel.  The  assumption  of 
vanishing  saturation  was  abandoned;  instead,  a  water  flux  was  set. 
This  allows  a  high  saturation  level  in  the  case  of  a  low  water  out¬ 
flow.  To  describe  the  saturation  at  the  interior  boundary  GDL  CL, 
a  continuous  saturation  was  assumed  contrary  to  the  suggestion  of 
a  continuous  capillary  pressure  [33,32].  This  results  in  a  saturation 
distribution  with  a  maximum  value  in  the  catalyst  layer. 

This  paper  presents  a  transient  PEM  fuel  cell  model  based  on  the 
assumptions  of  [39],  where  the  GDL,  CL  and  membrane  are  spatially 
resolved  in  1 D  with  an  agglomerate  approach  for  the  structure  of 
the  CL.  Unlike  the  model  of  Shah  et  al.  [38],  we  assume  a  discon¬ 
tinuity  of  the  saturation  due  to  the  continuous  capillary  pressure 
and  introduce  immobile  saturation  due  to  the  mixed  wettability 
of  the  GDL  structure.  The  anode  is  assumed  as  non-polarizable 
since  hydrogen  is  used  as  fuel,  therefore  the  anode  only  acts  as 
sophisticated  interface  for  the  water  management.  The  model  is 
validated  with  experimental  data  for  current  response,  impedance 
and  temperature  evolution  during  voltammetry  measurements. 


2.  Model  description 

The  present  model  is  developed  to  investigate  the  dynamic  tran¬ 
sients  of  a  PEMFC  in  detail.  For  this  purpose,  all  components  which 
are  responsible  for  a  characteristic  response  with  a  certain  time 
constant  are  considered  in  the  model.  These  are  namely  the  cath¬ 
ode  gas  diffusion  layer  (GDL)  that  suffers  flooding  in  the  case  of 
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Fig.  1.  Schematic  diagram  of  the  model  domains.  The  white  arrows  depict  the  domains  where  the  solving  variables  are  defined.  The  catalyst  layer  (CL)  consists  of  a  network 
of  agglomerates  where  the  oxygen  reduction  reaction  take  place.  When  water  is  present  in  the  CL,  it  is  assumed  that  the  water  covers  the  agglomerates  in  the  form  of  a  thin 
water  film  which  causes  an  additional  mass  transport  resistance  for  the  reactant. 


insufficient  water  removal,  the  cathode  catalyst  layer  (CCL)  where 
the  three-phase  boundary,  the  activation  overpotential  and  flood¬ 
ing  have  a  strong  influence  on  the  oxygen  reduction  reaction  (ORR) 
and  the  membrane  with  a  conductivity  that  is  highly  sensitive  to 
its  water  content  and  is  thus  responsible  for  ohmic  losses.  Due  to 
the  fast  kinetics  of  the  hydrogen  oxidation  and  the  high  hydrogen 
diffusivity,  the  polarization  losses  on  the  anode  side  are  neglected 
in  the  model.  We  treat  the  anode  as  a  spatially  reduced  interface  for 
the  boundary  conditions.  To  obtain  a  realistic  temperature  distribu¬ 
tion  within  the  membrane  electrode  assembly  (MEA),  the  latter  is 
sandwiched  between  two  model  domains  that  serve  for  the  bipolar 
plates  (BPs). 

Altogether,  the  model  consists  of  five  model  domains  (£2)  which 
are  schematically  depicted  in  Fig.  1.  All  layers  are  normalized  to 
a  thickness  of  one  for  both  the  numerical  applications  and  better 
visualization.  The  real  thickness  is  denoted  by  L Q . 

2  A.  Considering  one  agglomerate 


part  of  the  agglomerates  which  is  the  ionomer  phase.  For  mathe¬ 
matical  simplification,  it  is  assumed  that  the  water  covers  the  entire 
agglomerate  and  forms  a  thin  water  film  surrounding  the  agglom¬ 
erates  with  a  layer  thickness  d.  This  implies  an  additional  transport 
barrier  for  the  oxygen.  The  oxygen  first  has  to  dissolve  in  water, 
then  diffuse  through  the  water  film  to  reach  the  ionomer  of  the 
agglomerate  where  the  reaction  can  take  place  by  diffusing  to  the 
active  sites. 

The  diffusion  of  dissolved  oxygen  is  described  by  Fields  law: 


Jlo2{'\  =  -D'o2cdo2,s 


3Co2[r] 

9r  ’ 


(1) 


where  j’q2  is  the  oxygen  flux  in  the  medium  v  (v=water  agglomer¬ 
ate),  Dq2  is  the  effective  oxygen  diffusion  coefficient  in  the  medium 
v  and  Cq2  the  normalized  dissolved  oxygen  concentration  referred 
to  the  dissolved  oxygen  concentration  on  the  water-film  surface 
Cq2  s  (interface  void-water  film): 


Experimental  data  of  Ihonen  et  al.  [40]  suggest  an  agglomerated 
structure  in  the  catalyst  layer  with  a  bimodal  pore  size  distribu¬ 
tion.  Carbon  particles  supporting  the  catalyst  (10-20  nm  radius) 
agglomerate  to  form  clusters  (agglomerates),  which  are  filled  with 
ionomer  acting  as  a  binder.  The  resulting  ionomer-flooded  pores 
with  a  radius  of  about  3-10  nm  (called  primary  pores)  enable  proton 
conductivity  and  diffusion  of  dissolved  reactants  in  the  agglomer¬ 
ates.  The  larger  pores  between  neighboring  agglomerates  are  called 
secondary  pores  ( 10-50nm).  The  ORR  take  place  inside  the  agglom¬ 
erates  where  the  carbon-supported  catalyst  forms  a  three-phase 
boundary  together  with  the  ionomer  and  the  dissolved  oxygen.  It 
is  assumed  that  the  agglomerates  are  spherically  shaped  with  a 
mean  radius  Ra  and  homogeneously  distributed  in  the  CL.  The  con¬ 
tact  between  the  agglomerates  is  sufficient  to  ensure  the  free  flow 
of  protons  and  electrons.  The  tortuosity  of  the  flow  path  is  included 
in  the  calculation  by  a  Bruggeman  correction.  Before  the  oxygen 
reaches  the  catalyst  in  the  agglomerates,  it  has  to  dissolve  in  the 
ionomer  and  diffuse  to  the  active  sites.  In  the  presence  of  liquid 
water  in  the  secondary  pores  of  the  CL,  water  covers  the  hydrophilic 


(2) 


The  distance  from  the  center  of  the  agglomerate  is  denoted  by  the 
spherical  coordinate  r  in  brackets.  The  dissolved  oxygen  concentra¬ 
tion  on  the  water- film  surface  is  calculated  from  the  gaseous  oxygen 
concentration  c§2  s  above  the  water  film  in  the  secondary  pores  of 
the  CL  using  Flenry’s  law: 


lo2,s 


=  Hcgn 


02  ,s  ’ 


(3) 


where  H  is  known  as  the  Flenry  constant  [41  ].  The  effective  oxygen 
diffusion  coefficient  in  the  agglomerate  is  calculated  by  the 
Bruggeman  correlation: 


na  _  ^1.5  r\i 
u0  2~ea  u02 ’ 


(4) 


where  ea  is  the  fraction  of  primary  pores  filled  with  ionomer  and 
Dq2  is  the  diffusion  coefficient  of  dissolved  oxygen  in  the  ionomer. 
The  diffusion  coefficient  of  dissolved  oxygen  in  the  water  film  is 
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used  from  bulk  water  diffusion  at  60°  C  (in  m2  s-1)  according  to 
Bird,  Stewart  and  Lightfoot  [42]: 

D£2  =4.82  X  10“9.  (5) 

The  oxygen  mass  balance  is  given  by  Eq.  (6).  Oxygen  is  not  con¬ 
sumed  in  the  water  film,  whereas  the  oxygen  is  reduced  in  the 
agglomerate  by  the  ORR  described  by  Tafel  kinetics  [21,22,24],  a 
simplified  approach  of  the  Butler-Volmer  kinetics  justified  at  high 
overvoltages  (which  is  the  case  at  a  cell  voltage  of  900  mV  and 
below): 

0  if  Ra  <  r  <  Ra  +  d, 

-k0Xc^  sc^2  [r]eanFri/RT  if  0  <  r  <  Ra, 

(6) 


solving  variable  of  the  model)  can  be  determined  by  the  following 
calculation. 

Considering  a  unit  cell  with  the  volume  of  one  spherical  agglom¬ 
erate  (Vq  =  1 7 rR2)  and  the  surrounding  void  space  (secondary 
pores)  per  agglomerate  Vp ry,  the  porosity  of  the  CL  can  be  deter¬ 
mined  by 


cCL 


V^  +  Va 


From  this  it  follows: 


dry  V« 

”  ~  1  ' 


(13) 


(14) 


where  k0  is  the  reaction  rate,  a  the  transfer  coefficient,  n  the  number 
of  transfered  electrons,  R  the  gas  constant  and  T  the  local  temper¬ 
ature.  Following  Andreaus  et  al.  [43]  and  Vielstich  et  al.  [44],  both 
suggesting  a  dependence  of  the  transfer  resistance  and  reaction 
rate,  respectively,  on  the  ionomer  water  content  and  water  activity, 
respectively,  we  introduce  a  linear  dependency  of  the  ORR  rate  on 
the  ionomer  water  content  X  (definition  see  Section  2.3). 

The  following  boundary  conditions  are  used  to  calculate  the 
current  production  per  agglomerate. 


•  The  concentration  of  dissolved  oxygen  on  the  surface  of  the  water 
film  is  related  to  the  local  gaseous  oxygen  concentration  in  the 
secondary  pores  by  Henry’s  law  and  here  assumed  as  known.  The 
normalized  concentration  is  set  to  one: 


Cq2  [Ra  +  d]  —  1 . 


(7) 


•  A  continuous  oxygen  concentration  at  the  interface  water  film 
[R+]*> agglomerate  [R~]: 

cdo2K]  =  cdo2Kl  (8) 

•  A  continuous  oxygen  flux  at  the  interface  water  film 
[R+  ]  agglomerate  [R“  ]: 


Jo2K+]=jS2[Ra-]. 

•  A  vanishing  oxygen  flux  in  the  center  of  the  agglomerate: 


(9) 


(10) 


An  analytical  solution  for  the  oxygen  flux  jo2[Ra,y]  at  the 
agglomerate  surface  (r  =  Ra)  as  a  function  of  the  local  overpotential 
p[y],  water  film  thickness  d[y]  and  dissolved  oxygen  concentration 
on  the  surface  c ^  s[y]  can  be  found.  Using  Faraday’s  law,  the  current 
generation  of  one  agglomerate  can  be  calculated  by  integrating  the 
oxygen  flux  over  the  agglomerate  surface  95: 

jgenM  =  j)4Fjo2[Ra,y]dA 

as 

4jtRq4FCq2  s[y]  Dq2  (1  +  d[y])(-l  + 'I'coth['I']) 

"  d[y]  D«’  4/ cothl*]  +  djyj (D^  -  ’ 


where 


^  =  ^4 


/  k0XeoinF^/RI 


(12) 


is  known  as  the  Thiele  modulus  [45]. 

The  relationship  between  the  local  water-film  thickness  d[y]  and 
the  saturation  s[y]  of  the  secondary  pores  in  the  CL  (which  is  a 


The  local  saturation  of  the  secondary  pores  in  the  CL,  solved  in  Eq. 
(41 ),  is  expressed  by  the  ratio  of  the  volume  of  the  thin  water  film 
VWf  to  the  secondary  pore  volume  in  the  dry  state: 


s=  v 

v**’ 


(15) 


where 


Vwf=fj((Ra+d)3-R3a). 


(16) 


Inserting  Eqs.  (14)  and  (16)  into  Eq.  (15),  we  obtain  a  relationship 
between  the  water  film  thickness  (surrounding  the  agglomerates) 
and  the  saturation  of  the  CL  at  the  local  y-position: 


d[y]  =  Ra 


1) 


(i  +(s[y]-i) 


\-€CL 

1 


(17) 


Combining  Eqs.  (17)  and  (11 )  results  in  the  current  generation  of  an 
agglomerate  jgen[y]  as  a  function  of  the  local  overpotential,  gaseous 
oxygen  concentration  and  saturation.  The  analytical  expression  for 
Jgenly]  is  used  in  the  following  as  source/sink  terms  for  the  conti¬ 
nuity  equations  of  the  oxygen  concentration  and  the  overpotential 
in  the  homogeneous  model  of  the  cathode  CL  that  are  described  in 
the  next  sections. 


2.2.  Activation  and  ohmic  overpotential 

The  activation  overpotential  p  is  the  driving  force  for  the  ORR 
and  is  defined  as  the  deviation  of  the  difference  of  the  electronic 
potential  Oe  and  protonic  potential  from  the  theoretical  open 
circuit  potential  AO°  =  1.23  V: 

rj  =  AO°-(Oe-OP).  (18) 

The  overall  overpotential  of  a  fuel  cell  is  the  sum  of  the  activation 
overpotential  in  the  catalyst  layers  and  the  ohmic  overpotentials 
due  to  finite  conductivity  of  the  fuel  cell  components  and  contact 
resistances  at  the  interfaces  of  the  different  layers. 

Our  model  takes  into  account  the  activation  overpotential  in  the 
cathode  CL,  proton  migration  in  the  cathode  CL  and  in  the  mem¬ 
brane.  The  local  change  in  electric  potential  in  the  CL  is  negligible 
compared  to  the  change  of  the  protonic  potential  due  to  a  much 
higher  electric  conductivity  of  the  carbon  matrix  than  the  proton 
conductivity  of  the  ionomer  and  is  thus  assumed  to  be  constant 
through  the  CL.  The  anode  activation  overpotential  is  also  neglected 
due  to  the  fast  kinetics  of  the  hydrogen  oxidation  reaction.  The 
above  described  ohmic  overvoltages  of  the  electrons  in  the  GDLs 
and  BPs  as  well  as  the  contact  resistances  are  summarized  into  one 
model  domain,  more  precisely  into  the  anode  bipolar  plate  (ABP), 
and  modeled  with  a  constant  conductivity  orcontact • 
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Ohm’s  law  is  used  for  the  description  of  the  charge  flux: 


oQ  dp 

Jp’e  =  ~laty’ 


(19) 


where  the  superscript  Q  stands  for  the  model  domain  (CL,  mem¬ 
brane  and  ABP).  The  proton  conductivity  of  the  ionomer  phase 
depends  on  the  water  content  A  and  is  defined  according  to  Springer 
et  al.  [1]: 


f  eh5  (0.514  A  -  0.326)  e1268^/303-!/1)  if  £2  =  {LCL,Lmem} 

[^  G contact  if  £2  =  LABP . 

(20) 


For  the  ionomer  volume  fraction  in  the  membrane,  a  value  of  e,-  =  1 
is  used  and  in  the  catalyst  layer  =  cp,  whereas  a  Bruggeman 
correction  accounts  for  the  reduction  of  the  pathway  for  the  protons 
in  the  CL. 

The  charge  balance  equation  reads: 

%£=  IV  (-<?ORR  -Cdl^J  if  £2  =  La  (21 } 

dy  [o  if  n  =  {Lmem,LBP2}, 

where  Cdl  is  the  double  layer  capacity  and  q0RR  is  the  current  gener¬ 
ation  per  agglomerate  jaggi0  multiplied  by  the  agglomerate  density 
A  in  the  CL: 


QORR  =jgen  ■  A.  (22) 

The  agglomerate  density  is  determined  by  the  agglomerate  volume 
and  the  void  space: 


1  _  3  (1  —  CpL ) 

Va  +  V^~  4  nRl 


(23) 


2.3.  Water  content 


The  water  content  of  the  ionomer  A  is  defined  as  the  ratio  of  the 
number  of  water  molecules  to  the  number  of  charge  sites  (SO“ ) 
in  the  cathode  catalyst  layer  and  in  the  membrane.  The  flux  of 
dissolved  water  in  the  ionomer  is  determined  by  diffusion  and 
electro-osmotic  drag: 


e'-SpiD'  9 A  Oi ^ag  . 

EWLq  dy  + 

diffusion  electro-osmotic  drag 


(24) 


where  EW  and  pi  denote  the  equivalent  weight  and  the  density  of 
the  ionomer,  respectively.  The  diffusion  is  forced  by  the  gradient  of 
the  dissolved  water  concentration  in  terms  of  the  water  content.  It 
is  assumed  that  the  diffusion  coefficient  is  a  function  of  the  water 
content  according  to  [46]: 

D[  =  2.1  X  10“7  Ae“(2436/T).  (25) 

The  electro-osmotic  drag  coefficient  in  the  second  term  on  the 
right-hand  side  is  defined  as 


equilibrium  water  content  Aeq  (Eq.  (29))  and  the  desorption  in  the 
inverse  case  follows: 


Qad  ~ 


kads  Cy  Cv  (heq  —  A)(l— S)  if  A  <  Aeq, 

kdes  fijy  (1  -  S)  (Aeq  -  A)  (1  -  RH )  if  A  >  keq  A  RH  <  1 , 

(28) 


where  kads  and  kdes  are  the  adsorption  and  desorption  constants, 
respectively.  The  equilibrium  water  content  is  a  function  of  the 
water  activity  aw  given  by  Springer  et  al.  [1  ]: 

Aeq[aw]  =  0.043  +  17.81aw  -  39.85a^  +  36. 0a^,  (29) 


where  the  water  activity  is  equal  to  the  relative  humidity  RH,  cal¬ 
culated  by  the  vapor  concentration  (cv  •  cj)  (defined  in  section  2.6) 
and  the  saturation  pressure  psat[T]  for  a  given  temperature: 


aw  =  RH  = 


(cv  ■  c*)RT 
psat[T] 


(30) 


with  a  saturation  pressure  (in  Pa)  according  to  [1  ]: 


logio 


psat[r] 

101325 


-2.1794  +  0.02953  (T  -  273.15) 


-9.1837  x  10-5  (T  -  273. 15)2 
+ 1.4454  x  10-7(T-  273. 15)3.  (31) 


According  to  the  cluster-network  model  of  Weber  and  Newmann 
[35,36],  we  assume  a  maximum  water  uptake  of  the  ionomer  up 
to  A™ax  =  14  (given  by  Eq.  (29)  at  RH  =  100%)  if  the  membrane 
is  in  contact  with  water  vapor.  In  the  case  of  liquid  water  on 
the  membrane  surface,  ionomer  restructuring  due  to  expanding 
hydrophobic  channel  occurs,  resulting  in  a  maximum  water  con¬ 
tent  of  A™ax  =  22.  The  fraction  of  expanded  channels  and  thus  the 
water  content  is  a  strong  function  of  the  existing  capillary  pres¬ 
sure  p[  in  the  ionomer.  A  relationship  is  suggested  by  Weber  and 
Newmann  [36]: 

A[pc]  =  X™*  +  1(A.|"»-A™“) 

ln[-1.6(rw  cos[(9m]]  -  ln[10 5 plc] 

0.3  V2 


1  -  erf 


where  erf  is  the  error  function.  The  uptake  of  liquid  water  by  the 
ionomer  if  the  capillary  pressure  p^L  in  the  secondary  pores  of  the 
CL  (defined  in  Section  2.4)  is  higher  than  in  the  pore  network  of  the 
ionomer  plc  and  the  liquid  water  release  from  the  ionomer  in  the 
inverse  case  is  given  by 


/  MPc-P“)  if  Pc<P“> 
\kr(Pc-P“)  if  Pc  >  PcL- 


(33) 


It  is  assumed  that  the  ORR  produces  water  in  dissolved  form 
Qorr/ 2F.  A  water  accumulation  term  accounts  for  the  time  depen¬ 
dence.  The  source/sink  terms  are  confined  in  the  CL  and  vanish  in 
the  membrane. 


2.4.  Saturation 


_  2.5A 

ttdrag  -  -22_- 


(26) 


The  mass  conservation  equation  describes  the  change  in  the  dis¬ 
solved  water  flux: 


Oh 

dy 


lq 


—  Qur  + 


Qorr 
2  F 


Cj  Pi  dA\ 
EW  dt  J  ’ 


(27) 


coupled  with  phase  transition  and  water  generation.  The  adsorp¬ 
tion  of  water  vapor  into  the  ionomer  at  a  water  content  below  the 


The  saturation  s  is  defined  in  the  GDL  and  CL  as  the  volume  of 
liquid  water  divided  by  the  volume  of  secondary  pores  in  the  dry 
state.  A  high  saturation  reduces  the  flow  pathway  for  the  oxygen  in 
the  porous  media  that  finally  results  in  an  increased  mass  transport 
resistance.  The  transport  of  the  liquid  water  in  a  porous  media  is  not 
completely  understood  but  it  is  well  known  that  the  flow  in  (both 
hydrophilic  and  hydrophobic)  porous  media  occurs  always  from  the 
higher  saturation  region  towards  the  lower  saturation  region.  The 
gradient  of  the  capillary  pressure  is  mostly  stated  as  the  dominant 
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Fig.  2.  ESEM  image  of  a  Toray  TGP-H-090  GDL  shows  partly  hydrophilic  regions, 
leading  to  the  model  assumption  of  immobile  saturation. 


immobile 


Fig.  3.  A  scaled  Leverett  Function  is  applied  for  the  capillary  pressure-saturation 
relationship. 


force  for  the  water  transport.  We  use  Darcy’s  law  for  modeling  the 
liquid  water  transport  in  porous  media: 


Kabs  Krel  dpc 

LQ  /x  dy  ’ 


(34) 


where  I<^s  is  the  absolute  permeability  and  /x  is  the  water  vis¬ 
cosity.  The  relative  permeability  Krei  and  the  capillary  pressure  pc 
are  sensitive  functions  of  the  saturation.  There  is  no  consense  in 
the  description  of  the  functional  variation  of  the  capillary  pres¬ 
sure  with  saturation  pc[s].  The  most  common  approaches  are  the 
Leverett  J-function,  the  Van  Genuchten  model  [47]  and  the  Brooks- 
Corey  model  [48,49]  which  all  have  their  origin  in  hydrology 
(ground  water  simulation,  water/oil  flow  in  unconsolidated  rock, 
etc.).  Natarajan  et  al.  [27]  suggest  a  simple  exponential  approach, 
agreeing  with  observed  polarization  curves. 

Environmental  scanning  electron  microscope  (ESEM)  images  of 
a  Toray  TGP-H-090  clearly  show  a  mixed  wettability  of  the  fiber 
structure  with  hydrophobic  regions  and  hydrophilic  regions  (Fig.  2). 
This  leads  to  the  assumption  of  immobile  saturation  (sim)  in  our 
model  description.  Starting  from  a  completely  dry  GDL,  gener¬ 
ated  liquid  water  from  the  CL  first  fills  all  the  hydrophilic  parts 
of  the  adjacent  GDL  pores.  If  these  hydrophilic  pores  are  filled,  the 
local  saturation  exceeds  sim  and  the  water  moves  by  capillary  force 
towards  the  lower  saturation  region,  where  again  the  hydrophilic 
pores  get  filled  first.  If  all  hydrophilic  regions  of  the  GDL  are  filled 
with  liquid  water  then  a  mobile  phase  continuity  from  CL  to  the  gas 
channel  is  reached.  Below  the  threshold  of  sim,  water  occupying 
the  hydrophilic  pores  has  to  evaporate  for  saturation  to  decrease 
further.  For  the  expression  of  pc[s]  we  followed  the  most  widely 
employed  form,  the  Leverett  J-function,  but  considering  the  immo¬ 
bile  saturation,  in  the  modified  form  suggested  in  [33]: 


p f?  =  crw  cos[(9^] 


(35) 


whereby  the  saturation  of  the  empirical  J-function,  originally  pro¬ 
posed  by  Udell  [30],  is  replaced  by  the  reduced  saturation  s*: 


•S  i — > 


S  Sim  . 
1  —  Sim 


Ja[s]  =  1 .417  s,  -2.12  si  +  1.263  sj. 


(36) 

(37) 


The  difference  to  the  original  approach  can  be  seen  in  Fig.  3. 

We  assume  a  continuous  capillary  pressure  pc  at  the  inter¬ 
face  between  CL*>  GDL.  Due  to  the  different  pore  size  distribution, 


porosity  and  wetting  properties  (contact  angle)  of  the  two  media 
(CL/GDL),  a  continuous  capillary  pressure  results  in  a  discontinu¬ 
ous  saturation  distribution  across  the  interface  [33,32].  Using  the 
chain  rule  for  the  derivation  of  the  capillary  pressure  dpc/dy  = 
(d pc/ds){ds/dy),  Darcy’s  law  (Eq.  (34))  can  be  written  as  a  simple 
diffusion  equation  for  the  saturation: 


.  __D?ds 
Js~  L°dy' 


(38) 


where  the  diffusion  coefficient  itself  is  a  function  of  the  saturation: 


Kabs  Krel  dp £ 

5  /x  ds 


(39) 


with  a  relative  permeability  following  the  cube  of  saturation: 

Krel  =  (s-Simf.  (40) 


The  continuity  equation  yields: 


^p  djs 

vw  dy 


Ln 


+  dur  -  — 
Vw 


(41) 


where  the  interfacial  mass-transfer  rate  between  water  vapor  and 
liquid  water  q^c  is  assumed  to  be  proportional  to  the  saturation  in 
the  porous  media  Q  and  proportional  to  the  difference  between 
the  existing  relative  humidity  and  fully  humidified  gas: 

„  _{  keva^stfcKRH-  1)  if  RH  <  1, 

Hec  —  \  Vw  l^z) 

{  kcon(l  -s)e^c;c„(RH-  1)  if  RH  >  1 . 

For  relative  humidity  ( RH )  below  100%,  liquid  water  evaporates 
(upper  term)  while  for  supersaturated  gas,  water  vapor  condenses 
(lower  term).  The  evaporation  rate  keva  and  the  condensation  rate 
kcon  are  assumed  to  be  constant  [50]. 


2.5.  Oxygen  concentration 


The  oxygen  concentration  is  defined  in  the  CL  and  GDL.  We  use 
simple  Fick  diffusion  for  the  flux: 


*5, 

L  Q  dy  ’ 


(43) 


where  is  the  normalized  gaseous  oxygen  concentration  referred 
to  the  oxygen  concentration  at  the  boundary  to  the  gas  channel 
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c^* ,  calculated  by  the  ideal  gas  law  with  the  known  oxygen  partial 
pressure  po2  at  the  cell  inlet  and  the  channel  temperature  Tch : 


?*  _  Po2 
°2  RTch 


(44) 


The  values  of  the  GDL  porosity  published  in  fuel  cell  models  are 
prevalently  below  0.5  [20,25,13,51,19].  According  to  the  manu¬ 
facturer,  common  GDLs  have  a  porosity  higher  than  0.7  in  the 
uncompressed  state  [52].  The  discrepancy  between  these  two  val¬ 
ues  cannot  only  be  explained  by  the  high  compression  level  in  the 
assembled  cell.  The  answer  can  be  found  in  the  model  assumptions. 
A  single-phase  model  (1, 2  or  3D)  does  not  account  for  liquid  water 
occupying  a  part  of  the  void  space  which  results  in  a  lowered  effec¬ 
tive  porosity  [53,25].  A  two-phase  ID  model  needs  a  low  porosity 
as  a  model  parameter  to  account  for  in-plane  inhomogeneities  not 
only  over  the  entire  MEA  but  also  on  the  scale  of  a  few  millimeters 
as  a  result  of  shadowing  in  the  rib^>channel  geometry.  A  simple  dif¬ 
fusion  simulation  in  2D  shows  that  the  shadowing  by  the  rib  results 
in  minimized  apparent  porosity  for  a  reduced  ID  model. 

The  effective  oxygen  diffusivity  accounts  therefore  not  only  for 
the  porosity  (with  a  realistic  value)  and  the  local  saturation  but  also 
for  the  spatial  reduction  to  a  1 D  model  with  the  geometry  factor  E: 


D^°=(S6«(1-s))1-5D*2.  (45) 

The  temperature  and  pressure  dependence  of  the  free  oxygen  dif¬ 
fusivity  follows  the  Chapman-Enskog  formula  [42]: 


Taking  the  reaction  into  account,  the  mass  balance  yields: 


^4 

dy 


<fO 


(47) 


where  adsorption/desorption  qad ,  evaporation/condensation  qec 
and  a  water  vapor  accumulation  term  are  included. 


2.7.  Temperature 


The  heat  flux  is  described  by  conduction  in  the  solid  media  and 
by  convection  in  the  form  of  liquid  water,  defined  in  all  five  layers. 
Heat  convection  by  the  gas  flux  is  neglected.  The  resulting  equation 
for  the  heat  flux  is 

h  =  ~lMy+  ’  (53) 

S  ^  convection  by  liq .  water 

conduction 

where  kq  is  the  thermal  conductivity  of  the  medium  £2  and  Q  = 
Pi  q  is  the  product  of  the  density  of  the  species  i  and  the  specific 
heat  capacity.  The  energy  equation  reads: 


djr 

dy 


T  £2 


nQ  rn  dr 
dt 


(54) 


where  the  source  term  qff  stands  for  Joule  heating,  latent  heat  or 
reaction  and  activation  loss,  depending  on  the  model  domain  £2. 

In  the  cathode  bipolar  plate,  Joule  heating  by  the  cell  current  icen 
due  to  the  finite  electrical  conductivity  a contact  is  assumed: 


qCBP=LCBP 


cell 


4  O contact 


(55) 


where  the  total  electrical  and  contact  resistance  (see  Section  2.2) 
is  split  into  components  for  the  cathode  BP  (25%),  the  cathode  GDL 
(25%)  and  the  anode  BP  (50%,  including  the  anode  GDL). 

Joule  heating  and  latent  heat  caused  by  the  phase  transition  in 
the  GDL  yields: 


2. 6.  Wa ter  vap or  concert tra don 


The  governing  equation  for  the  vapor  flux  jv  in  the  porous  media 
is  given  by 


D f'aCj  dcv 
dy  ’ 


(48) 


where  cv  is  the  normalized  vapor  concentration  referred  to  the  inlet 
vapor  concentration  c*.  The  latter  is  calculated  by  the  saturation 
pressure  psat  (Eq.  (31))  for  a  given  humidification  temperature  T^p 
of  the  cathode  humidifier  applying  the  ideal  gas  law: 


c  — 


Psat[TcDP] 

RT 


(49) 


The  effective  diffusion  coefficient  accounts  for  the  porosity  of 
the  porous  media  by  the  Bruggeman  correlation  and  the  clogging  of 
liquid  water,  which  hinders  vapor  diffusion.  The  factor  E  accounts 
for  the  shadowing  effects  of  the  rib  in  a  real  fuel  cell  (see  Section 
2.5): 


neff& 

Uy 


(S€"(l -S))1'5^ 


(50) 


where  Df  is  the  temperature-dependent  diffusivity  in  a  non-porous 
condition  according  to  Um  and  Wang  [54]: 


D*  =  7.35  x  1(T5 


(51) 


The  mass  conservation  equation  can  be  expressed  as 


_tQ 

dy 


^  Qad  Qec 


1  —  S)Cy 


dCp\ 

■dfj- 


(52) 


% 


V  4  (7 contact 


+  hgi  qec 


(56) 


In  the  CL,  the  reaction  heat  of  the  ORR  is  the  dominant  factor,  but 
also  latent  and  ohmic  heat  is  included: 


tf=La 


( AST  +  Frj)  +  hgiqec 


reaction  heat 


latent  heat ) 


ohmic  heat 


(57) 


Proton  migration  in  the  membrane  produces  heat  due  to  the  finite 
protonic  conductivity: 

Cm  =  ~Jp  (58) 

The  heat  production  in  the  anode  bipolar  plate  corresponds  to  that 
of  the  cathode  BP: 

qABP=LABP  'Jell  (59) 

2  CT contact 

The  heat  capacity  CQ  of  certain  domains  £2  are  listed  below,  where 
the  newly  introduced  subindices  (Ti  and  C)  stand  for  titanium  and 
carbon: 

CBP  :  CCBP  =  CTi, 

GDL  :  CCDL  =  e^DLsCw  +  e™t(i  _ s)cair  +  (1  -  e£Di)Cc. 

CL  :  CCL  =  ecpLsCw  +  e“(l  -  s)Cair  +  (1  -  ecpL)Cc,  (60) 

Mem  :  Cmem  =  Cmem 
ABP  :  CABP  =  CTi. 
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membrane  CL  GDL  channel 


Fig.  4.  The  schematic  diagram  of  a  spatially  resolved  anode  illustrates  the  vapor  distribution  in  the  cases  of  desorption  (dashed)  and  adsorption  (solid).  The  reduced  anode 
interface  accounts  for  ad-/desorption  and  water  vapor  diffusion,  implemented  as  a  lumped  boundary  condition. 


2.8.  Boundary  conditions 

We  refer  to  Fig.  1  for  the  discussion  of  the  boundary  conditions 
where  the  solving  variables  and  their  area  of  validity  are  indicated.  If 
nothing  else  is  stated,  we  assume  continuity  in  the  solving  variable 
and  continuous  flow  at  the  interior  boundaries. 


As  for  oxygen,  the  membrane  is  also  impermeable  for  the  water 
vapor: 


dcv[  1] 
dy 


(66) 


2.8 A.  Overpotential  r\ 

Since  the  protons  are  not  allowed  to  penetrate  into  the  GDL,  their 
flux  at  this  interface  (y  =  0)  is  taken  to  be  zero: 

«-0,  (6,) 

dy 

At  the  boundary  of  the  anode  bipolar  plate  (y  =  3 ),  the  overpotential 
is  set  to  a  value  dependent  on  the  simulated  cell  voltage  U: 


f?[3]  =  1.23V-  U-rja,  (62) 

where  r]a  is  the  anodic  overpotential  which  is  not  calculated  in 
the  model  but  can  be  extracted  from  the  experimental  data  for 
reference  electrodes  [55]. 


2.8.2.  Oxygen  concentration  c§2 

The  average  oxygen  concentration  at  the  interface  GDL  **  flow 
channel  depends  strongly  on  the  stoichiometry  and  flow  rate  due 
to  oxygen  consumption  along  the  channel.  To  take  the  air  flow  rate 
Vc  into  account,  a  Cauchy-type  boundary  condition  is  chosen: 


f02[- i]  =  vc^2cg;a-cg2).  (63) 

It  gives  a  relationship  between  the  air  flow  rate,  the  deviation  of  the 
normalized  oxygen  concentration  from  100%  and  a  specific  oxygen 
transfer  coefficient  Qc0^ ,  chosen  large  enough  that  Cq2  never  falls 
below  90%,  even  for  a  high  load  (which  is  realistic  for  small  test  fuel 
cells). 

We  assume  that  the  membrane  is  an  impermeable  barrier  for 
gaseous  oxygen  which  implies  zero  flux  at  y  =  1 : 


dy 


=  0. 


(64) 


2.8.3.  Water  vapor  concentration  cv 

The  flux  boundary  condition  for  the  water  vapor  concentration 
is  defined  similarly  to  the  oxygen  concentration  described  above: 

jv[-\\  =  VcQcvc*v{\-cv\  (65) 


2.8.4.  Water  content  X 

Since  dissolved  water,  in  terms  of  X,  has  no  transport  mechanism 
in  the  gas  diffusion  layer,  its  flux  is  taken  to  be  zero  at  the  interface 
between  CL  and  GDL: 


dX[0] 

~dy~ 


(67) 


The  boundary  condition  for  the  water  content  on  the  anode  side  is 
extremely  complex  and  needs  a  detailed  explanation. 

In  consideration  of  the  fact  that  on  the  anode  side,  neither  the 
catalyst  layer  nor  the  gas  diffusion  layer  is  spatially  resolved,  the 
relevant  transport  mechanism  and  phase  transitions  therein  have 
to  be  reduced  to  an  interface  with  a  sophisticated  boundary  condi¬ 
tion.  For  better  understanding  of  the  layer  reduction  to  a  simplified 
interface,  an  illustration  of  the  local  variables  is  given  in  Fig.  4.  In 
order  to  allow  dehydration  of  the  ionomer  on  the  anode  side  we  do 
not  fix  the  water  content  A,  to  a  certain  value  but  use  a  Cauchy-type 
boundary  condition  as  follows. 

We  calculate  the  equilibrium  water  vapor  partial  pressure  at  the 
membrane  surface  paymem  for  a  given  water  content  at  the  anode 
interface  (A.[2])  by  means  of  the  inverse  function  XXq[aw]  of  Eq. 
(29)  and  the  saturation  pressure  psat[T ]: 


pOymem  =psat[T]x-l[aw].  (68) 

The  vapor  pressure  can  be  converted  by  means  of  the  ideal  gas  law 
to  a  vapor  concentration: 


a, mem 
c  v 


a, mem 
Pv 

RT 


(69) 


Neglecting  vapor  diffusion  in  the  catalyst  layer,  we  use  the  CL 
thickness  for  a  simple  adsorption/desorption  region,  where  a  linear 
dependence  of  the  adsorption/desorption  process  on  the  difference 
between  the  equilibrium  vapor  concentration  Cy,mem  at  the  mem¬ 
brane  surface  and  the  existing  vapor  concentration  c„,GDL  between 
CL  and  GDL  is  assumed.  The  vapor  flux  at  the  boundary  CL  **  GDL 
can  be  written  as 


where  Qcv  is  the  specific  water  vapor  transfer  coefficient  at  the 
interface  GDL  **  channel. 


za  tCL  j.a  r  ra,  GDL  .a, mem  s 
)ad  =  L  Kad^Cv  ~cv  h 


(70) 
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where  lCL  is  the  anode  CL  thickness  and  kaad  the  adsorp¬ 
tion/desorption  rate,  expressed  as 


/5xl0-5kdes  if  cav’mem  >  c°-CDL, 
\5xl0  ~Hads  if  c°;mem  <  c°-GDL. 


(71) 


The  vapor  diffusion  in  the  GDL,  forced  by  the  vapor  concentration 
gradient,  is  written  as 


2.8.6.  Temperature  T 

With  regard  to  our  water-cooled  fuel  cell  (which  provided  the 
measured  data  we  use  for  the  model  validation),  we  use  Cauchy 
boundary  conditions  that  consider  a  linear  dependence  of  the  heat 
removal  on  the  difference  between  the  temperature  at  the  outer 
side  of  the  BPs  and  the  cooling- water  temperature: 


M- 2]  =  jr[3]  =  £2t(T  -  Tcooiant ), 


(79) 


{„a,cha  _a,i 

\a  _  nefT&  { Cv  ~cv 

Jdiff  ~  "  iGDl 


where  T2t  is  the  heat  transfer  coefficient  dependent  on  the  coolant 
(72)  flow  rate. 


where  is  the  effective  vapor  diffusion  coefficient  and  Cy,cha  is 
the  vapor  concentration  in  the  anode  channel. 

The  condition  of  continuity  of  both  fluxes,  provided  that  there 
is  no  phase  change  in  the  GDL,  results  in 


fad  =  Jdiff  = 


-Df,  (c"' 


,mem  _  a,cha 
Lv 


a)lCLka 


ad 


Dy  +  LCL  ka  ,  LGDL 


(73) 


As  a  boundary  condition  for  the  water  content  on  the  anode  side, 
we  couple  the  dissolved  water  flux  with  the  derived  vapor  flux: 


h 


-D?  (c° 


mem  _  a,cha 
—  Lv 


)Lcl  k“ 


ad 


Df,  +  LCL  kan.  LGDL 


C2°w(\-u)en-u], 

vg 


(74) 


where  the  last  term  on  the  right-hand  side  accounts  for  the  case 
of  high  water  content  {X  >  14)  when  the  ionomer  releases  liquid 
water,  carried  out  by  the  hydrogen  gas  stream.  A  high  value  of  fi*,  = 
600  mr2  is  chosen  for  fast  water  release  and  6[x]  is  the  Heavyside 
function  (equal  to  1  for  x  >  0  and  0  elsewhere). 


2.9.  Impedance 

The  membrane  impedance  (Z)  is  calculated  by  an  integral  over 
the  inverse  membrane  conductivity  (Eq.  (20))  along  the  membrane 
thickness: 

C2  T  mem 

z  =  l  ^dy-  (80) 

2.10.  Numerical  details 

The  governing  equations  are  solved  using  COMSOL  Multiphysics 
™,  a  commercial  software  package  based  on  finite  element  meth¬ 
ods.  A  Direct  Linear  System  Solver  (UMFPACK)  and  quadratic 
Lagrange  polynomials  as  test  functions  are  used.  The  solver  is 
allowed  to  take  free  time  steps.  The  five  model  domains  are  dis¬ 
cretized  with  a  non-uniform  grid  of  265  elements  whereas  the  CL 
and  the  regions  near  the  interfaces  between  the  domains  are  refined 
with  smaller  elements. 


2.8.5.  Saturations 

The  accumulation  and  removal  of  liquid  water  on  top  of  the 
GDL  (towards  the  channel)  is  a  complicated  process.  Water  forms 
droplets  which  are  blown  out  by  the  convective  force  of  the  gas 
stream.  The  greater  the  flow  velocity  in  the  channel,  the  more  effi¬ 
cient  the  removal  of  liquid  water.  Such  a  sophisticated  dynamic 
approach  can  hardly  be  reproduced  in  a  ID  model.  As  a  strongly 
simplified  approach,  we  implemented  an  outflow  condition  for  the 
saturation  as  a  quadratic  function  of  the  difference  between  the 
saturation  itself  and  the  immobile  saturation  sim : 

js[-l]  =  /  -^nUs~sim)2  lf  S  >  Sjm ,  (75) 

lo  ifs<sim, 

where  the  flow  rate  has  a  linear  influence  on  the  liquid  water  flux. 

The  interior  boundary  of  the  saturation  differs  from  the  other 
solving  variables  by  having  a  discontinuity  in  its  distribution.  Con¬ 
tinuous  capillary  pressure  is  assumed  at  the  interface  GDL  CL: 

p?*[0]=p?[0],  (76) 

which  is,  in  our  opinion,  the  physically  realistic  boundary  condi¬ 
tion.  This  results  in  a  saturation  discontinuity  at  the  interface  (in 
the  case  of  different  structural  properties  and  wettabilities  such  as 
porosities  and  contact  angles)  but  with  a  continuous  liquid  water 
flux: 


J?D1[0]  =  jsa[0]. 


(77) 


A  zero-flux  boundary  for  liquid  water  is  assumed  at  the  interface 
CL  membrane: 


asm 

dy 


=  0. 


(78) 


3.  Experimental 

The  experiments  were  made  with  a  small  test  fuel  cell  (geo¬ 
metric  area  of  1  cm2 )  in  a  test  bench  that  are  both  described  in 
detail  in  Gerteisen  [55].  Note  that,  due  to  the  reference  electrode 
configuration  used,  we  are  able  to  separate  the  anode  overvoltage 
from  the  remaining  cell  losses,  which  is  important  for  the  valida¬ 
tion  of  a  fuel  cell  model  where  the  anode  losses  are  neglected.  An 
untreated  gas  diffusion  layer  of  type  TORAY®TGP-H-090  was  uti¬ 
lized.  A  GORE  ™  PRIME  A®  Series  5510  MEA  with  Pt  catalyst  (loading 
anode/cathode:  0.4  mg  cm-2)  and  a  membrane  thickness  of  35  p,m 
was  used.  For  analyzing  water  transport,  flooding  effects  and  mem¬ 
brane  dehydration,  dynamic  measurements  were  performed  that 
give  a  good  insight  into  the  processes  that  occur.  We  used  voltam¬ 
metry  and  chronoamperometry  experiments.  To  prove  the  validity 
of  the  model,  we  investigated  the  current  response  of  the  cell 
for  three  different  humidification  conditions  of  the  inlet  gases, 
denoted  as  case(l):  airdly/H2diy,  case  (2):  airdry/H2hum  and  case  (3): 
airhum/^ 2hum  -  Detailed  operating  conditions  such  as  the  dew  point 
temperature  of  the  inlet  gases  and  flow  rates  are  summarized  in 
Table  1. 

Voltammetry  experiments  were  conducted  in  the  voltage 
range  between  900  mV  and  60  mV  with  a  scan  rate  of  10mVs-1. 


Table  1 

Operating  conditions  in  the  voltammetry  experiments 
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current  density  /  Acm'2 


Several  sweeps  were  applied  sequentially  to  reach  and  ensure 
reproducibility  of  the  dynamic  characteristic  of  the  I-U  curve.  The 
cell  impedance  at  10  kHz  was  recorded  during  all  experiments.  The 
temperature  was  measured  with  a  thermocouple  in  the  anode  plate 
1  mm  above  the  flow  field  . 

Chronoamperometry  is  a  measurement  technique  in  which  the 
cell  voltage  is  switched  between  two  values,  and  the  resulting  cell 
current  response  is  monitored  as  a  function  of  time.  The  charac¬ 
teristic  of  the  current  response  and  cell  impedance,  triggered  by  a 
voltage  change  from 

•  step  (1):  600-300  mV, 

•  step  (2):  700-400  mV, 

•  step  (3):  800-500  mV  and  vice  versa,  was  investigated. 

4.  Results  and  discussion 

4.1.  Voltammetry  experiments 


Fig.  5.  A  typical  voltammetry  measurement  that  shows  MEA  dehydration  at  high 
cell  voltage  and  flooding  in  the  limiting  current  density  region  which  results  in  two 
hysteresis  loop  like  a  “eight”.  By  means  of  reference  electrode  measurements  the 
overall  cell  loss  can  be  separated  into  cathode,  anode  and  ohmic  losses. 


Voltammetry  experiments  were  used  to  validate  the  developed 
model  with  measured  data.  Therefore,  we  briefly  discuss  the  exper¬ 
imental  results  to  establish  an  understanding  of  the  measured 
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Fig.  6.  Voltammetry  experiments  with  a  scan  rate  of  10  mVs-1  were  conducted  to 
validate  the  developed  model,  (a)  Voltammetry  experiments  conducted  under  dif¬ 
ferent  humidification  conditions  show  strong  hysteresis  effects,  (b)  Measured  cell 
impedances  show  a  large  dynamic  range  in  dehydration  and  humidification  of  the 
ionomer. 


Fig.  7.  The  simulated  polarization  curves  and  impedance  characteristics  show  a 
similar  influence  on  the  humidification  conditions  as  the  measurement,  (a)  The 
simulated  current-voltage  characteristic  show  strong  hysteresis  effects  with  about 
the  same  limiting  current  density  as  the  measurement,  (b)  The  increase  of  the 
impedance  lies  in  the  same  range  as  the  experiment. 
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dynamic  effects.  As  an  example,  Fig.  5  shows  the  cell  voltage,  sepa¬ 
rated  into  cathode,  anode  and  ohmic  overvoltage  during  a  dynamic 
voltammetry  experiment  (forward  and  backward  sweep)  for  dry 
operating  conditions.  It  can  be  seen  that  the  anode  overvoltage  has 
minor  impact  on  the  overall  cell  loss  but  is  not  completely  negligi¬ 
ble,  which  has  to  be  in  mind  for  the  validation.  Fig.  6(a)  shows  the 
polarization  curves  of  voltammetry  experiments  for  three  different 
humidification  conditions  (cases  (l)-(3)).  The  corresponding  cell 
impedance  is  depicted  in  Fig.  6(b).  The  forward  (900-60  mV)  and 
the  backward  sweeps  (60-900  mV)  are  labeled  by  arrows.  For  com¬ 
parison  with  the  simulation  results,  the  cell  voltage  is  corrected  by 
the  anode  losses  and  denoted  by  the  superscript  *.  The  polarization 
characteristic  is  strongly  influenced  by  the  humidification  condi¬ 
tions  in  such  a  way  that  under  completely  dry  conditions  (case  (1 )), 
the  highest  limiting  current  density  is  reached  at  the  expense  of 
ionomer  dehydration  at  low  current  density.  Fig.  6(b)  shows  for  case 
( 1 )  that  the  impedance  increases  throughout  a  backward  sweep  and 
further  increases  in  the  forward  sweep  up  to  a  maximum  value  of 
approximately  230  m  £2cm2 ;  thereafter  the  water  production  is  suf¬ 
ficient  to  humidify  the  ionomer  itself.  The  increasing  impedance  in 
the  low  current  density  range  is  responsible  for  the  hysteresis  effect 
in  the  polarization  curve  in  that  region.  The  observed  lower  limiting 
current  density  in  the  case  of  higher  humidification  (cases  (2)  and 
(3)),  leads  to  a  decline  from  1.6Acm_2to  1.3  Acm_2and  0.9  A  cm-2, 
respectively,  indicating  higher  oxygen  transport  limitations  due  to 
increased  flooding  phenomenon.  The  saturation  level  during  the 
voltammetry  experiment  is  strongly  time-dependent  and  causes 
hysteresis  loops  in  the  high  current  density  region  in  all  three  cases. 
The  dehydration  of  the  ionomer  is  considerably  reduced  in  case  (2), 
compared  to  the  completely  dry  conditions  in  case  (1 ),  and  almost 
prevented  for  case  (3). 

4.2.  Model  validation 

The  simulation  results  for  the  three  different  sets  of  operating 
conditions  are  shown  in  Fig.  7.  The  values  of  the  model  parameters 
used  for  the  simulation  are  listed  in  the  parameter  list  in  Section 
6  and  are  also  used  in  simulations  that  are  discussed  in  the  next 
sections. 

Under  all  three  different  operating  conditions,  the  model 
predicts  the  same  current-voltage  characteristics  as  the  measure¬ 
ments,  as  shown  in  Fig.  7(a).  In  all  three  cases  the  simulations  show 
hysteresis  effects  in  the  high  current  density  region  with  about 
the  same  limiting  current  density.  In  case  (1),  the  pronounced  hys¬ 
teresis  loop  in  the  low  current  density  region  is  observed  too.  The 
time-dependent  behavior  of  the  simulated  impedance  (Fig.  7(b)) 
is  in  good  agreement  with  the  measurement.  In  case  (1)  with 
completely  dry  conditions,  the  simulated  impedance  shows  also 
a  strong  variation  with  the  same  maximum  value  as  the  experi¬ 
ment.  In  case  (2),  the  dehydration  of  the  ionomer  is  less  pronounced 
and  the  impedance  in  case  (3)  is  nearly  constant,  as  found  in  the 
measured  data. 

A  comparison  of  the  temperature  change  throughout  a  sweep 
experiment  is  made  in  Fig.  8.  On  the  left-hand  side,  the  measured 
temperature  is  shown  for  the  three  different  sets  of  humidifica¬ 
tion  conditions.  It  is  evident  that  the  cell  cooling  is  not  sufficient 
to  keep  the  temperature  constant  during  a  potential  sweep.  The 
high  current  density  in  case  (1)  is  responsible  for  a  temperature 
rise  of  3.5  K  from  313.5  K  to  317  K.  For  comparison,  we  analyze 
the  simulated  temperature  at  the  interior  boundary  between  the 
anode  plate  and  the  anode  interface  T[2],  depicted  in  Fig.  8(b). 
The  simulated  temperature  shows  the  same  qualitative  charac¬ 
teristic  as  the  measurement.  The  temperature  variation  falls  in 
the  same  range  of  3.5 1<  for  case  (1),  3K  for  case  (2)  and  2I<  for 
case  (3). 


current  density  /  Acrrf2 


current  density  /  Acm'2 

Fig.  8.  A  qualitative  agreement  between  the  measured  temperature  in  the  anode 
plate  and  the  simulated  temperature  is  achieved,  (a)  The  temperature  of  the  cell 
increases  up  to  2-3.5  K  depending  on  the  load,  (b)  The  non-isothermal  model  pre¬ 
dicts  qualitatively  the  same  dynamics  for  the  cell  temperature. 

4.3.  Analysis  of  the  solving  variables 

In  this  subsection,  we  discuss  the  origin  of  the  hysteresis  effects 
in  detail  on  the  basis  of  the  profiles  of  certain  solving  variables.  The 
following  plots  show  snapshots  of  some  profiles  either  at  0.7  V  or 
0.4  V  during  a  simulated  voltammetry  experiment. 

For  the  analysis  of  the  hysteresis  loop  in  the  low  current  den¬ 
sity  region,  which  is  most  pronounced  in  case  (1),  we  refer  to  the 
water  content  X  at  0.7  V.  In  Fig.  9(a),  a  large  difference  in  the  water 
content  between  forward  and  backward  sweeps  is  present  in  cases 
(1)  and  (2).  The  low  water  content  in  the  forward  mode  results 
in  a  high  ohmic  drop  in  the  membrane  which  finally  results  in 
a  minimized  activation  overpotential  in  the  cathode  CL  for  the 
ORR  (Fig.  9(b)).  In  addition,  the  lower  water  content  reduces  the 
area  of  the  three-phase  boundary  zone  and  therefore  decelerates 
the  reaction  kinetics  (see  Eq.  (6)).  The  fact  that  oxygen  transport 
limitation  is  not  responsible  for  the  lowered  current  generation, 
which  might  be  possible  in  the  case  of  different  saturation  levels, 
is  shown  in  Fig.  10.  In  the  forward  mode,  more  oxygen  is  avail¬ 
able  for  the  ORR  than  in  the  backward  mode.  Fig.  11(b)  shows 
that  the  relative  humidity  in  the  CL  remains  at  nearly  100%  except 
for  forward  sweeps  of  cases  (1)  and  (2)  where  dry  cathode  inlet 
gases  are  used.  The  difference  of  saturation  between  forward  and 
backward  sweeps  at  0.4  V  is  depicted  in  Fig.  11(a).  The  saturation 
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voltage:  0.7  V 


voltage:  0.7  V 


Fig.  9.  A  snapshot  of  the  profiles  of  the  water  content  and  the  overpotential  between 
forward  and  backward  sweeps  at  0.7  V  is  given,  (a)  A  comparison  of  the  water  con¬ 
tent  profile  shows  a  strong  dehydration  of  the  ionomer  in  cases  (1)  and  (2).  (b) 
The  steep  ohmic  drop  in  the  membrane  in  case  (1)  leads  to  a  minimized  activation 
overpotential  in  the  CL,  resulting  in  low  current  generation. 

level  in  a  backward  sweep  is  always  above  its  level  in  the  forward 
mode.  The  largest  difference  is  observed  in  case  (1 ),  where  the  sat¬ 
uration  differs  up  to  20%  in  the  GDL.  The  saturation  level  of  the 
GDL  in  our  simulations  is  uncommonly  high  compared  to  previ¬ 
ous  publications  [51,3,33,8,10,11].  The  reason  for  it  is  based  firstly 
on  the  assumption  of  an  existing  immobile  saturation  of  sim  =  0.2, 
and  secondly  in  the  choice  of  Cauchy  boundary  conditions  at  the 
interface  GDL  channel.  The  most  common  assumption  for  the 
boundary  condition  is  vanishing  saturation  at  the  interface  towards 
the  channel.  The  high  capillary  force  in  the  presence  of  a  large  sat¬ 
uration  gradient,  predicted  by  the  Leverett  J-function,  results  in 
a  fast  liquid  water  outflow.  This  implies  a  small  curvature  of  the 
saturation  profile  and  finally  a  low  saturation  level  in  the  GDL  in 
the  case  of  a  Dirichlet  boundary  condition:  s  =  0  at  the  interface 
GDL  channel.  Wang  et  al.  [56]  predicts  a  maximum  saturation 
level  of  only  about  6.3%  at  a  current  density  of  about  1.4  A  cm-2  in 
the  case  of  the  commonly  made  assumptions.  Such  a  low  saturation 
level  would  not  lead  to  mass  transport  limitations  at  a  current  den¬ 
sity  of  about  1.0  A  cm-2  (using  realistic  values  of  the  GDL  porosity), 
that  are  clearly  observed  in  our  experiments. 


voltage:  0.7  V 


Fig.  10.  A  comparison  of  the  profiles  of  the  normalized  oxygen  concentration  (0.7  V) 
during  a  forward  and  a  backward  sweep.  Higher  transport  limitations  are  observed 
in  the  backward  sweep. 

To  the  best  of  our  knowledge,  Shah  et  al.  [39]  are  among  the  few 
using  a  Cauchy  boundary  condition,  leading  to  a  saturation  of  up  to 
0.6-0.8.  Unlike  in  our  model,  they  assume  continuity  of  the  satura¬ 
tion  at  the  interface  CL  GDL  and  therefore  reached  the  maximum 
saturation  level  in  the  CL. 

Fig.  12(a)  shows  the  time-dependent  system  response  in  terms 
of  current  density  and  saturation  at  the  interface  GDL  ^channel 
for  case  (1).  Both  the  maximum  of  current  density  (A<pi )  and  sat¬ 
uration  (A  cp2)  display  a  phase  shift  relative  to  the  minimum  of  the 
applied  cell  voltage.  The  cumulative  saturation  during  the  forward 
mode  causes  increasing  mass  transport  limitation,  which  becomes 
the  dominant  loss  mechanism  in  the  high  current  density  region. 
This  in  turn  leads  to  a  current  density  maximum  which  is  consider¬ 
ably  earlier  than  the  voltage  minimum  (about  25  s).  The  saturation 
increases  beyond  this  voltage  minimum  up  to  a  phase  shift  of  about 
Acp2  =  90°  (45  s).  The  evolution  of  the  ionomer  water  content  at  the 
interface  GDL  CL  is  analyzed  in  Fig.  12(b).  The  maximum  of  the 
water  content  shows  only  a  small  phase-shift  relative  to  the  voltage 
(A  cp3),  whereas  the  minimum  of  the  water  content  shows  a  strong 
delay  with  respect  to  the  voltage  maximum  (A <p4). 

The  evolution  of  the  water  content  profile  for  case  (2)  is  visual¬ 
ized  in  Fig.  13.  In  the  region  of  maximum  current  density,  labeled  as 
(1 ),  a  water  content  of  nearly  17  is  reached  in  the  CL.  A  strong  gra¬ 
dient  of  the  water  content  is  observed  in  the  y-direction  due  to  the 
strong  electro-osmotic  drag  at  high  current  densities.  The  decreas¬ 
ing  electro-osmotic  drag  during  the  backward  sweep  and  the  back 
diffusion  leads  to  a  balancing  of  the  water  content  in  the  membrane 
(region  (2)).  In  region  (3),  the  introduced  water  from  the  humidi¬ 
fied  hydrogen  stream  at  the  anode  and  the  water  generation  from 
the  ORR  in  the  cathode  CL  is  insufficient  to  prevent  dehydration  of 
the  ionomer.  The  water  content  falls  to  a  value  of  7.1. 

4.4.  Chronoamperometry 

A  comparison  between  measured  and  simulated  chronoamper¬ 
ometry  experiments  is  made  in  Figs.  14  and  15.  The  fuel  cell  was 
operated  with  dry  air  and  humidified  hydrogen  (Tgp  =  311  K)  with 
the  same  flow  rate  as  in  the  voltammetry  experiments  (see  Table  1 ). 
The  simulated  voltage  steps  are  corrected  by  the  measured  anode 
losses.  Since  the  anode  overvoltage  shows  small  changes  in  time 
after  a  voltage  step,  the  mean  value  is  used.  Overshoots  in  the 
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voltage:  0.4  V 


voltage:  07  V 


Fig.  11.  Comparison  of  the  profiles  of  the  saturation  (0.4  V)  and  the  relative  humidity 
(0.7  V)  during  a  forward  and  a  backward  sweep,  (a)  Continuity  in  the  capillary  pres¬ 
sure  at  the  interface  results  in  a  discontinuity  in  the  saturation.  A  higher  saturation 
is  found  in  the  backward  mode,  (b)  The  relative  humidity  (RH)  for  case  (3)  remains 
at  nearly  100%  during  the  entire  sweep. 
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Fig.  12.  The  time-dependent  characteristic  of  the  applied  cell  voltage  together  with 
the  system  response  in  terms  of  cell  current,  saturation  (s[ — 1  ])  and  water  content 
(A.[0])  is  shown  for  case  (1 ).  (a)  The  maximum  current  density  is  reached  about  20  s 
before  the  lowest  cell  voltage  is  reached  (A<pi ).  The  maximum  saturation  (interface 
channel  **  GDL)  has  a  delay  ( A(p2 )  of  about  42  s  relative  to  the  minimum  cell  voltage, 
(b)  The  characteristic  of  the  water  content  (interface  GDL  CL)  follows  that  of  the 
current  density  in  principle.  However,  the  decrease  of  the  water  content  is  a  little 
delayed  (A (p4)  with  respect  to  the  cell  current  due  to  the  liquid  water  reservoir  in 
the  catalyst  layer. 


measured  current  density  are  recorded  on  switching  from  high  to 
low  cell  voltages  (Fig.  14(a)).  The  characteristic  of  the  overshoots 
depends  on  the  voltage  level  before  and  after  the  step  change.  A 
change  from  0.6  to  0.3  V  (step  (1))  causes  a  peak  current  density 
of  1.3  Acm_2that  drops  off  to  a  value  of  0.9  A  cm-2  within  the  first 
20  s.  The  peak  current  density  of  step  (2)  reaches  the  same  value 
but  shows  a  less  rapid  decay  than  step  (1).  The  peak  in  case  (3)  is 
not  very  pronounced,  but  shows  an  apparent  delay  relating  to  the 
step  change,  that  is  not  observed  for  step  (1)  and  (2)  to  this  degree. 

Switching  from  the  low  to  the  high  cell  voltage,  the  characteristic 
of  the  current  response  differs  for  all  three  steps.  Step  (1)  shows 
an  increase  in  the  current  density  immediately  after  the  voltage 
change  and  reaches  equilibrium  within  100  s.  Step  (2)  shows  a  short 
increase  or  plateau  before  the  current  density  starts  to  decrease.  The 
current  density  in  step  (3)  decreases  continuously  to  a  steady-state 
value. 

The  simulation  shows  a  qualitatively  similar  characteristic 
(Fig.  14(b)).  Flowever,  the  time  constants  of  the  response  differ  from 
the  measured  ones.  The  simulation  of  steps  (1)  and  (2)  predicts  a 


distinct  overshoot  in  the  current  density  with  the  same  peak  value 
of  1.5  A  cm-2.  A  sharp  peak  is  observed  for  step  (1),  but  step  (2) 
shows  a  much  slower  decay  than  the  measured  data.  Step  (3)  shows 
a  linear  decrease  of  the  current  without  a  clear  peak.  A  delay  of  the 
maximum  current  relative  to  the  step  change  is  predicted  by  the 
model,  in  agreement  with  the  measurement.  The  simulated  cur¬ 
rent  response  to  the  potential  step  from  low  to  high  cell  voltage 
is  in  agreement  with  the  measurement.  An  increase  of  the  current 
density  for  step  (1),  an  increase  in  the  current  density  followed  by 
decrease  for  step  (2)  as  well  as  the  fast  decay  to  a  steady-state  value 
for  step  (3)  are  predicted. 

The  corresponding  cell  impedance  was  recorded  during  the 
chronoamperometry  measurements  and  is  depicted  in  Fig.  15(a). 
Almost  no  change  in  the  impedance  was  measured  for  step  (1 ).  The 
ionomer  starts  to  dehydrate  when  the  fuel  cell  is  driven  at  a  cell  volt¬ 
age  of  700  mV  (step  (2)).  Stepping  back  to  400  mV,  the  impedance 
drops  very  rapidly  to  a  low  impedance  level,  similar  to  the  value  of 
step  (1 ).  At  a  cell  voltage  of  800  mV,  the  impedance  rises  very  fast 
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Fig.  13.  The  3D  diagram  visualizes  the  evolution  of  the  water  content  profile  during 
two  cycles  of  the  sweep  experiment  for  case  (2).  The  water  content  is  dominated  by 
electro-osmotic  drag  and  water  generation  (region  1 ),  back  diffusion  (region  2)  and 
dehydration  (region  3). 


Fig.  14.  Comparison  between  measured  and  simulated  current  response  for 
chronoamperometry  experiments,  (a)  An  overshot  is  observed  in  all  three  cases  on 
switching  to  the  low  cell  voltage.  The  response  differs  on  returning  to  the  higher 
voltage,  depending  on  the  voltage  level,  (b)  A  qualitatively  good  agreement  with  the 
experiments  was  achieved  in  the  simulations,  but  with  a  longer  time  constant. 


time/s 

Fig.  15.  Comparison  between  measured  and  simulated  impedance  responses  during 
chronoamperometry  experiments,  (a)  An  increase  in  the  impedance  at  cell  voltages 
of  700  mV  and  above  shows  that  the  self-humidification  by  product  water  is  insuf¬ 
ficient.  (b)  The  simulation  predicts  dehydration  in  the  voltage  range  of  700  mV  and 
above,  but  the  simulated  impedance  is  not  as  high  as  the  measured  impedance. 

to  a  value  of  0.35  £2 cm2  and  drops  again  very  fast  on  switching  to 
higher  current  density  and  overvoltage. 

The  simulation  of  the  impedance  response  describes  the 
measured  impedance  qualitatively.  While  switching  between 
600-300  mV,  dehydration  of  the  ionomer  can  not  be  observed.  In 
the  range  between  700-400  mV,  the  ionomer  starts  to  dehydrate 
at  a  cell  voltage  of  700  mV.  Strong  dehydration  of  the  ionomer  at 
a  cell  voltage  of  800  mV  for  step  (3)  is  predicted  by  the  model,  but 
is  obviously  not  as  pronounced  as  in  the  measured  data.  A  reason 
for  this  could  be  the  shrinking  of  the  ionomer  in  the  case  of  high 
dehydration,  resulting  in  an  increase  of  the  contact  resistance  that 
is  not  implemented  in  the  model. 

5.  Conclusion 

In  PEM  fuel  cells,  water  is  present  as  vapor,  in  liquid  form 
and  dissolved  in  the  ionomer.  The  effects  of  water  in  all  three 
phases,  basically  pore  flooding  and  ionomer  dehydration,  are  inves¬ 
tigated  with  a  newly  developed  ID  transient  model.  The  model 
accounts  for  the  loss  mechanisms  in  the  cathode  GDL,  cathode 
CL  and  membrane.  The  electrode  structure  is  modeled  as  a  net¬ 
work  of  spherical  agglomerates.  An  excess  of  liquid  water  in 
the  catalyst  layer  leads  to  coverage  of  these  agglomerates,  form¬ 
ing  a  water  film  which  limits  oxygen  transport.  Liquid  water  is 
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modeled  as  saturation  in  the  void  space.  Its  transport  properties 
depend  strongly  on  the  capillary  pressure-saturation  relationship 
in  porous  media,  which  in  turn  is  a  function  of  the  wettability 
and  pore  structure.  Based  on  ESEM  images,  immobile  saturation 
is  introduced  due  to  the  observed  partly  hydrophilic  regions  in 
the  GDL.  We  assume  a  continuous  capillary  pressure  at  the  inter¬ 
face  CL  GDL,  resulting  in  a  discontinuous  saturation  distribution. 
Finite  phase  transition  rates  between  the  ionomer  and  pore  space, 
namely  adsorption/desorption  and  liquid  water  uptake/release, 
enable  membrane  dehydration  in  the  case  of  dry  operating  con¬ 
ditions,  resulting  in  an  increased  ohmic  resistance. 

The  simulation  results  for  dynamic  current-voltage  character¬ 
istics  are  in  excellent  agreement  with  our  measured  data.  The 
measured  hysteresis  loop  in  the  limiting  current  density  region 
is  reproduced  by  the  model  and  explained  by  pore  flooding.  The 
measured  cell  impedance  during  a  voltage  sweep,  which  is  a  mea¬ 
sure  for  the  time-dependent  water  content,  is  in  agreement  with 
the  simulation  over  a  wide  range  of  operating  conditions  and  cur¬ 
rent  densities.  Simulations  of  chronoamperometry  experiments  are 


performed  with  the  validated  parameter  set.  The  model  captures 
qualitatively  the  current  and  impedance  responses  of  measured 
chronoamperometry  data. 

Our  model  forms  a  suitable  basis  for  identifying  the  dominant 
loss  mechanisms  at  different  operating  points.  Based  on  the  simula¬ 
tion  results,  suggestions  for  a  better  choice  of  fuel  cell  components 
can  be  made,  guidance  for  their  physical  properties  like  pore  struc¬ 
ture,  wettability  or  thickness  can  be  obtained  or  simply  optimal 
operating  conditions  can  be  found.  The  latter  leads  us  to  the  topic  of 
fuel  cell  control.  A  good  understanding  of  the  response  on  dynamic 
load  changes  is  required  to  control  a  fuel  cell  within  a  stable  oper¬ 
ating  point.  This  can  be  realized  by  implementing  our  model  in  a 
simplified  form  in  a  model-based  control  algorithm. 

Parameter  studies  and  a  sensitivity  analysis  would  be  the  next 
steps  in  ongoing  work.  An  upgrade  to  a  multi-dimensional  model 
is  intended  for  the  future. 

6.  Nomenclature  and  parameter  list 


Symbol 

Description 

Value  /  Eq. 

Unit 

Reference 

Solving  variables 

o2 

Normalized  gaseous  oxygen  concentration 

- 

cv 

Normalized  vapor  concentration 

- 

s 

Saturation 

- 

T 

Local  temperature 

I< 

ri 

Overpotential 

V 

X 

Water  content 

- 

Physical  constants 

F 

Faraday  constant 

96484 

Cmol-1 

R 

Gas  constant 

8.314 

JK”1  mol”1 

Structural  values 

Lcbp 

Thickness  of  cathode  bipolar  plate 

3  x  10-3 

m 

Lgdl 

GDL  thickness 

280  x 10~6 

m 

Lgl 

CL  thickness 

10  x 10-6 

m 

jjnem 

Membrane  thickness 

35  x  10-6 

m 

Labp 

Thickness  of  anode  bipolar  plate 

3  x  10-3 

m 

Ra 

Mean  agglomerate  radius 

0.2  x  10-6 

m 

[57] 

£a 

Volume  fraction  of  primary  pores  in  an  agglomerate 

0.4 

- 

[26] 

eCL 

p 

Volume  fraction  of  secondary  pores  in  CL 

0.25 

- 

[40,22] 

€CL 

Volume  fraction  of  ionomer  in  CL 

0.3 

- 

Jgdl 

p 

Volume  fraction  of  open  pores  in  GDL 

0.7 

- 

[58] 

A 

Agglomerate  density 

Eq.  (23) 

m-3 

3 

Geometry  factor 

0.55 

- 

Calculated 

Physical  properties,  local  variables  and  boundary  conditions 

aw 

Water  activity 

Eq.  (30) 

- 

cd 

02,S 

Dissolved  oxygen  concentration  at  agglomerate 

Eq.  (3) 

mol  ITT3 

surface 

cd 

LOo 

Normalized  dissolved  oxygen  concentration  in 

Eq.  (2) 

- 

u2 

agglomerate 

£ 

Gaseous  oxygen  concentration  at  inlet 

Eq.  (44) 

mol  ITT3 

e* 

Vapor  concentration  at  inlet 

Eq.  (49) 

mol  ITT3 

Cdl 

Double  layer  capacity 

0 

Fm-3 

CTi  =  PTi  cTi 

Heat  capacity  of  titanium 

9.4  x  106 

J  ITT3  K”1 

[59] 

CGDL  =  pGDL  CGDL 

Heat  capacity  of  GDL 

1.61  x  106 

J  ITT3  K-1 

[38] 

CCL  =  Pa  CCL 

Heat  capacity  of  CL 

1.61  x  106 

J  ITT3  K-1 

[38] 

Cmem  _  pmem  cmem 

Heat  capacity  of  membrane 

2.18  x  106 

J  ITT3  K-1 

[38] 

Eg  —  Pg  cg 

Heat  capacity  of  gas 

1  x  103 

J  m-3  K_1 

[38] 

Gw  —  Pw  Cw 

Heat  capacity  of  water 

4.187  x  106 

J  m-3  K_1 

[38] 

d 

Water  film  thickness 

Eq.  (17) 

m 

Do2 

Oxygen  diffusion  coefficient  in  ionomer 

5  x  10-11 

m2  s_1 

[22] 

Do2 

Oxygen  diffusion  coefficient  in  water 

2.1  x  10-9 

m2  s_1 

Oxygen  diffusion  coefficient  in  gas  phase 

Eq.  (46) 

m2  s_1 

[42] 

peff,£2 

O? 

Effective  oxygen  diffusion  coefficient  in  gas  phase 

Eq.  (45) 

m2  s_1 

[42] 

D? 

Vapor  diffusion  coefficient  in  gas  phase 

Eq.  (51) 

m2  s_1 

nCj ff.fl 

Effective  vapor  diffusion  coefficient  in  gas  phase 

Eq.  (51) 

m2  s_1 

Df 

Water  diffusion  coefficient  in  porous  medium 

Eq.  (39) 

m2  s_1 

Di 

Water  diffusion  coefficient  in  ionomer 

Eq.  (25) 

m2  s_1 

[46] 

EW 

Equivalent  weight 

1.1 

kg 

[20] 

hgi 

Heat  of  vaporisation/condensation 

40.7  x  103 

J  mol-1 

[59] 

H 

Henry  constant 

0.0254 

- 

[41] 

180 
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Symbol 

Description 

Value  /  Eq. 

Unit 

Reference 

4 

Gaseous  oxygen  flux 

Eq.  (43) 

mol  m-2  s-1 

iv 

Vapor  flux 

Eq.  (48) 

mol  m~2  s-1 

is 

Liquid  water  interstitial  velocity 

Eq.  (34) 

ms-1 

)T 

Heat  flux 

Eq.  (53) 

Jrrr2  s_1 

Jp,e 

Charge  flux 

Eq.  (19) 

Am-2 

ik 

Dissolved  water  flux 

Eq.  (24) 

mol  m-2  s-1 

Jgen 

Current  generation  per  agglomerate 

Eq.  (11) 

ka , 

Adsorption/desorption  rate  on  anode 

Eq.  (71) 

s-1 

Assumed 

^ ads 

Adsorption  rate 

80 

s-1 

Assumed 

kdes 

Desorption  rate 

50 

s-1 

Assumed 

keva 

Evaporation  rate 

100 

ms  kg-1 

[50] 

kcon 

Condensation  rate 

100 

s-1 

[50] 

ku 

Water  uptake  rate 

0.1 

mol  s  kg-1  m~2 

Assumed 

kr 

Water  release  rate 

0.1 

mol  s  kg-1  rrr2 

Assumed 

ko 

Reaction  rate 

0.001 

s-1 

Assumed 

KGDL 

Absolute  permeability  in  GDL 

0° 

Cj 

X 

o 

N) 

m2 

[39] 

l<ci 

Absolute  permeability  in  CL 

1  X  10-13 

m2 

[39] 

Krel 

Relative  permeability 

Eq.  (40) 

n 

Number  of  transfered  electrons 

2 

- 

p? 

Capillary  pressure  in  domain  Q 

Pa 

psat 

Saturation  pressure 

Eq.  (31) 

Pa 

[1] 

Qorr 

Source/sink  term  due  to  ORR 

Eq.  (22) 

Am-3 

Qad 

Source/sink  term  due  to  adsorption/desorption 

Eq.  (28) 

mol  ITT3  s-1 

Qur 

Source/sink  term  due  to  liquid  water  uptake/release 

Eq.  (33) 

mol  ITT3  s-1 

Qec 

Source/sink  term  due  to  evaporation/condensation 

Eq.  (42) 

mol  m~3  s-1 

Qh 

Source/sink  term  of  the  heat 

Eqs.  (55)-(59) 

J  ITT3  s-1 

RH 

Relative  humidity 

Eq.  (30) 

- 

sGDL/CL 

Immobile  saturation 

0.2/0 

- 

Assumed 

a 

Transfer  coefficient 

0.45 

- 

01 drag 

Electro-osmotic  drag  coefficient 

Eq.  (26) 

- 

AS 

Enthalpy  change 

162.2 

J  mol”1  K-1 

[38] 

®m 

Contact  angle  in  ionomer  channels 

90.02 

a 

[36] 

Contact  angle  in  porous  media  (GDL/CL) 

105/95 

° 

«GDL 

Thermal  conductivity  of  GDL 

1.67 

Wm1  K-1 

[38] 

“CL 

Thermal  conductivity  of  CL 

0.67 

Wm1  K-1 

[38] 

Kmem 

Thermal  conductivity  of  membrane 

0.67 

Wm”1  K-1 

[38] 

KTi 

Thermal  conductivity  of  titanium 

21.9 

Wm”1  K-1 

[59] 

heq 

Equilibrated  water  content 

Eq.  (29) 

- 

[1] 

[X 

Liquid  water  viscosity 

0.001 

kgm-1  s_1 

[2] 

Vg 

Molar  volume  of  ideal  gas 

2.2414  x  10-2 

m3  mol-1 

Vw 

Molar  volume  of  liquid  water 

1.8015  x  10~5 

m3mol_1 

Pi 

Ionomer  density 

1980 

kg  m-3 

[20,18] 

crw 

Surface  tension  of  water 

0.07 

[2] 

G contact 

Unified  conductivity  of  all  electronic  conductors 

600/650 

Sm_1 

Assumed 

normed  to  BP2 

Protonic  conductivity  in  domain  £2 

Eq.  (20) 

Sm_1 

[1] 

°Co2 

Oxygen  transfer  coefficient  cathode 

2.4  x  104 

ITT2 

Assumed 

ai 

Vapor  transfer  coefficient  cathode 

2.4  x  104 

ITT2 

Assumed 

Qc 

Saturation  transfer  coefficient  cathode 

3  x  102 

ITT2 

Assumed 

Water  transfer  coefficient  anode 

1  x  10-5 

ITT2 

Assumed 

P2T 

Heat  transfer  coefficient 

1150 

J  K_1  m-2  s_1 

Assumed 

Operating  conditions 

hell 

Current  density 

A  cm-2 

Po2 

Oxygen  partial  pressure 

0.21 

atm 

P 

Total  pressure 

1 

atm 

T coolant 

Coolant  temperature 

313 

K 

' Td 

DP 

Anode  dew  point  temperature 

Various 

K 

Tc 

DP 

Cathode  dew  point  temperature 

Various 

K 

Va 

Anode  gas  flow  rate 

50(8.33  x  10-7) 

ml  min-1  (m3  s_1) 

Vc 

Cathode  gas  flow  rate 

100(16.66  x  10-7) 

ml  min-1  (m3  s_1) 

U 

Voltage 

V 
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